npj Systems Biology and Applications

Systems biology aims to develop predictive models that connect molecular mechanisms to cellular behavior. Genome-scale metabolic models are among the most widely used frameworks for integrating stoichiometric, thermodynamic, and omics-derived information to predict feasible metabolic phenotypes. However, cellular metabolism operates on timescales governed by enzyme kinetics and by the relationship between metabolic fluxes and metabolite pool sizes. In steady-state metabolic models, this relationship can be expressed in terms of metabolite turnover rates, defined as flux-to-pool-size ratios that quantify how rapidly metabolite pools are renewed. As a result, physiologically consistent steady-state solutions should not only satisfy mass-balance and thermodynamic constraints but also exhibit turnover rates consistent with enzyme-mediated cellular dynamics. Current constraint-based approaches can admit many steady-state flux-concentration states that do not account for turnover rates, resulting in phenotypes incompatible with realistic metabolic dynamics, even when multiple types of data are imposed. Here, we present METEOR-K, an optimization framework that links steady-state metabolic fluxes to metabolite concentrations via turnover rate constraints to identify dynamically plausible flux-concentration reference states. Because these constraints reshape the feasible solution space, we also introduce turnover-rate-aware sampling strategies to efficiently explore the resulting feasible region. We applied METEOR-K to models of increasing scope and scale, including a reduced glycolysis pathway, anaerobic E. coli, and near-genome-scale ovarian cancer models. METEOR-K narrowed the admissible steady-state solution space, reduced uncertainty in feasible flux-concentration states, and improved local dynamic behavior. In nonlinear ODE simulations of bioreactor cultivation and drug-response scenarios, METEOR-K-derived states produced intracellular response times compatible with growth-supporting metabolic operation and perturbation recovery. Overall, these results establish metabolite turnover rates as scalable biophysical constraints that improve the physiological consistency of steady-state metabolic modeling. Because turnover rates encode flux-to-pool-size timescale constraints, METEOR-K moves part of physiological-consistency assessment upstream of kinetic parameterization, yielding better-suited flux-concentration reference states for kinetic modeling and dynamic prediction.

Cancer cells rewire metabolism to support proliferation. Intriguingly, divergent metabolic choices are made to attain this common goal. Identifying the unique metabolic requirements for a specific cell has profound implications for cancer biology and precision medicine. Genome-scale metabolic models (GEMs) have emerged as powerful tools to systematically characterize, understand, and predict metabolism of cells and tissues. Despite being comprehensive, the current GEMs remain limited in their predictive power. Here, we present a new GEM of human cells, in silico Human Metabolic Essentiality (iHME), that significantly improves the prediction of metabolic dependencies at a reduced computational cost. Wse rationally downsized, curated, and corrected previous models to remove unsupported metabolic redundancies, which led to a slim model containing 4,377 reactions, 3,241 metabolites, and 1,825 genes. When used to reconstruct metabolic networks of 1,103 cancer cell lines, iHME recalled on average 84.6% of experimental essential genes, which is two-fold increase over previous models. Cholesterol biosynthesis was revealed to be the most reliably predicted pathway with alternative dependencies. Finally, we applied the model to reconstruct individualized networks and predict essential gene profiles for 8,384 patient tumor samples. Glucose transporter SLC2A1 (GLUT1) was identified as a context-specific dependency for head and neck cancers and ovarian cancer. Likewise, CDP-diacylglycerol synthase CDS2 was identified for skin cancer. Overall, iHME is a new genome-scale model for prediction of metabolic dependency at higher accuracy and computational efficiency.

Optimizing dose and schedule remains a central challenge in oncology drug development, particularly for immunotherapies where fixed dosing regimens often fail to account for patient specific heterogeneity in tumor-immune dynamics. Here, we present a hybrid quantitative systems pharmacology-reinforcement learning-Monte Carlo Tree Search (QSP-RL-MCTS) framework for personalized immunotherapy dosing that formulates dose selection as a sequential decision-making problem. The approach integrates a mechanistic QSP model of prostate cancer immunotherapy, transcriptomics informed virtual patient populations and data driven AI system comprising reinforcement learning and Monte Carlo tree search. Reinforcement learning is used to learn adaptive generalized dosing policies that optimize treatment outcomes across the population, while Monte Carlo Tree Search provides forward-looking evaluation of RL predicted dosing trajectories to refine patient-specific decisions. On benchmarking against fixed dosing regimens of ipilimumab, the remission rate of the proposed model (95.2%) was comparable to the highest fixed dosing regimen of 10 mg/kg per dose while the median total dose (72 mg/kg) of the proposed model designed regimen was comparable to the lowest fixed dosing regimen of 3 mg/kg per dose. The model is generalizable across different dosing protocols and can be extended to predict optimal dose under different therapeutic scenarios. Analysis of the learned dosing trajectories enables stratification of patients into distinct response groups and identifies drug activity rate as the dominant determinant of long-term treatment outcome. These results demonstrate how mechanistically guided artificial intelligence can transform population-level dose optimization into patient-specific, biologically interpretable treatment strategies for precision immuno-oncology.

Mechanistic ordinary differential equation models are widely used in systems biology to represent biochemical networks, population dynamics, cell-state transitions, and other biological processes; however, their predictive value depends critically on accurate parameter estimation from noisy and often sparse experimental data. In this tutorial, we present the Weak-form Estimation of Nonlinear Dynamics (WENDy) method as a forward-solver-free approach that reformulates parameter estimation as a covariance-corrected weak-form regression problem by integrating the model equations against compactly supported test functions. We present the background on the methodology through the lens of the familiar logistic equation, and we demonstrate applications of the method on real experimental data through two systems biology examples: a glycolytic oscillator with relatively dense time-course data and a sparse epithelial-mesenchymal cellstate transition model with multiple experimental replicates. Ultimately, using WENDy, we estimate interpretable biological parameters with uncertainty for systems with noisy and sometimes sparse available experimental data.

AO_SCPLOWBSTRACTC_SCPLOWCell fate decisions are driven by gene regulatory networks (GRNs). While the mutually inhibitory toggle switch effectively models binary fate decisions, fully connected inhibitory networks with more than two nodes fail to capture multi-fate decisions due to the low prevalence of "single high states", where only a single master regulator is highly expressed. The goal of this study is to find network structures that support all single high states. We find that the only network that attains the highest possible prevalence of all single high states within the set of monotone Boolean (MB) models is completely disconnected. Since biological networks typically require connectivity, we investigate network structures that support equipotency, where all single high states have equal prevalence within MB models. Finally, we characterize the networks that support multistability between all single high states, finding that it is possible only in networks in which each node either has self-activations or is inhibited by every other network node. Our findings provide a theoretical framework for understanding the network design principles that can support simultaneous differentiation into multiple distinct cell types.

Neuroblastoma (NB), a pediatric cancer arising from disrupted sympathetic neuron differentiation, exhibits marked heterogeneity and limited therapeutic options. To better understand its molecular circuitry dynamics, we applied CardamomOT, a novel Gene Regulatory Network (GRN) inference framework, to single-cell RNA-seq data from patient-derived tumoroids. This approach models gene regulation via piecewise deterministic Markov processes, capturing transcriptional bursting and protein-mediated feedback, overcoming limitations of RNA velocity (e.g., gene independence and lack of biological time). We identified a continuous chromaffin-to-sympathoblast differentiation trajectory along which we selected 85 dynamically relevant genes enriched in cell cycle and DNA replication functions. Notably, 9 genes overlapped with those driving normal sympathoadrenal differentiation, underscoring tumor-normal tissue similarity. The inferred 85-genes network reproduced quite well experimental gene expression patterns in silico, and allowed to predict protein-level dynamics. Furthermore, it allowed to predict the effect of perturbations (both knock-out and overexpression) of hub genes (e.g., tcf4 and PLK1). We show that those perturbations significantly altered cell fate proportions in silico, with tcf4 KO increasing chromaffin-like cells and reducing proliferative late sympathoblasts. Predictions regarding tcf4 were tested using drug inhibition as a proxy for the gene KO. Using the BET inhibitor JQ1 indeed induced profound effect on the transcriptomic identity of our tumoroids. All of the 50 predicted tcf4 target genes were found to be significantly altered by JQ1 treatment. Finally cell fate proportions were also altered ex vivo closely resembling the predicted output. Our work therefore demonstrates that NB tumoroids retain a dynamic, differentiation-like architecture amenable to GRN modeling. Predicted druggable targets offer testable therapeutic avenues, including repurposing BET inhibitors or PLK1 inhibitors, potentially in combination.

Multiomic datasets contain complex nonlinear relationships that are often missed by conventional analysis methods. Topological data analysis (TDA) can detect some such patterns (i.e., loops), and here we extend previous TDA approaches with MOLA (MultiOmic Loop Analysis), a framework for multiomic loop visualization, filtering, association with sample-level characteristics and functional characterization. In an Influenza A virus dataset, MOLA better captured anticipated characteristic-dependent epigenetic alterations compared to standard approaches. In a breast cancer cohort, higher loop participation in five genes was associated with increased hazard of progression. MOLA provides an interpretable framework for discovering biologically meaningful multiomic patterns and potential biomarkers.

Evaluating the potential applications of a medicine is a fundamental challenge in drug development. There is a lack of standardized, decision-oriented benchmarks that test whether computational models can generalize therapeutic hypotheses across diseases in ways that reflect real-world pharmaceutical investment decision making. To address this gap, we introduce two complementary resources: the Indication Expansion Investment Decision Network (IxIDN) and the Orphanet Rare Disease Ontology Negative-network (ORDON). IxIDN is a clinical-trial-derived positive benchmark constructed by projecting drug-disease associations from pharmaceutical clinical trials into a disease-disease network; each edge connects disease pairs that have entered clinical trials for the same drug, thereby capturing cases when concrete indication-expansion decisions have been made. The current release contains 574 rare diseases and 5,336 edges. In contrast, ORDON serves as a stringent, biology-aware negative benchmark derived from the authoritative Orphanet Rare Disease Ontology. It identifies maximally distant disease pairs according to curated hierarchical structure and genetics-linked inheritance patterns, providing 793 rare diseases and 5,000 edges that represent high-separation negative candidates across therapeutic areas. Together, IxIDN and ORDON enable rigorous cross-evidence generalization from clinical trials to disease ontology, testing for Disease-Disease Association Learning (DDAL), a core task for mechanism-centered drug repurposing and indication expansion. All data are publicly available with detailed metadata, enabling reproducible evaluation of models on transparent, decision-relevant benchmarks.

Computational modeling provides a powerful framework for in silico exploration of anti-cancer therapeutic targets and tumor response mechanisms. Oncogenic signaling pathways play a central role in tumor behavior and represent promising targets for personalized combination therapies. However, these pathways are complex, and although logic-based models are well suited for representing signaling dynamics, they are often constrained by model-specific data requirements, limited scalability, and time-consuming manual curation. Here, we introduce Functional Integration of Contextualized Omics for Unraveling regulatory dynamicS (FICUS), a framework that integrates omics-driven network contextualization with dynamic Boolean and logic-ODE modeling. FICUS enables automated, data-driven protein network inference and patient stratification, allowing shared signaling mechanisms to be identified across patient subgroups while preserving patient-specific dynamic responses. We applied FICUS to the SU2C-MARK lung cancer cohort and the The Cancer Genome Atlas kidney cancer cohort, demonstrating its utility for post-hoc analyses and downstream interrogation of dynamic tumor models. Overall, our results highlight the flexibility of FICUS in capturing heterogeneous signaling mechanisms across patient subgroups, addressing a key challenge in precision oncology.

We present DynaVelo, a generative neural ordinary differential equation model that learns the joint dynamics of gene expression and transcription factor (TF) motif activities in evolving cell systems using single-cell multiome with joint gene expression and chromatin accesibility readout. DynaVelo leverages partial RNA velocity information together with single-cell TF motif accessibility data to improve the modeling of cell state dynamics and identification of TF drivers. We show that DynaVelo recovers the complex and bifurcating in vivo dynamics of wildtype murine germinal center (GC) B cells and reveals how these cell dynamics change under loss-of-function mutations in epigenetic regulators Arid1a and Ctcf. DynaVelo resolves how TF motif activities evolve along latent time trajectories using analysis of training cells or through generated trajectories from the model. In silico perturbation analysis further enables DynaVelo to infer dynamic and cell-state-specific gene regulatory networks (GRNs), recovering many known TF-to-gene edges in the wildtype GC GRN and predicting those that are disrupted in mutants. Finally, in silico gene and TF perturbations allow both the prediction of cell dynamics under loss-of-function genetic mutations and the identification of TF perturbations to rescue loss-of-function dynamic and immunological phenotypes. This analysis predicted that Ctcf knockout would rescue Arid1a loss-of-function phenotype in the GC reaction and nominated Bcl6 and Stat3 as additional TFs whose knockout would rescue Arid1a loss. We validated these predictions in vivo using double heterozygous mutant mice, confirming rescue of the Arid1a dark zone phenotype in all cases and quantitatively assessing model predictions using multiome in Arid1aHet;CtcfHet double heterozygous mice. DynaVelo therefore provides a powerful new deep learning framework for modeling and perturbing dynamic cell systems by harnessing single-cell multiome data sets.

Spatial proteomic imaging technologies enable the simultaneous assessment of immune cell abundance and spatial organization within the tumor microenvironment. Spatial clustering is commonly summarized using measures such as Ripleys K or nearest-neighbor G-functions at a fixed radius. However, these approaches depend on scale selection and may obscure biologically relevant patterns occurring across spatial ranges. We propose a functional data analysis (FDA) framework to model spatial clustering trajectories derived across a continuum of radii. Functional principal component analysis (FPCA) was used to summarize dominant modes of spatial variation, and resulting scores were incorporated into Cox proportional hazards models as both main effects and interaction with immune cell abundance. The approach was applied to multiplex immunofluorescence data from five ovarian cancer studies, comprising 773 high-grade ovarian serous tumors. Analyses focused on CD3+ and CD8+ T cell populations within the tumor compartment of the tissue, adjusting for age at diagnosis and cancer stage, with study-specific estimates combined using random-effects meta-analysis. Higher abundance of both T cells and CD8+ T cells was consistently associated with improved overall survival. Beyond abundance, spatial features captured by the leading functional principal component were independently associated with survival, particularly for CD8+ T cells. Interaction models further showed that the prognostic effect of immune infiltration depended on spatial clustering, with tumors characterized by high abundance and low spatial clustering exhibiting the most favorable outcomes. These findings indicate that spatial organization provides complementary prognostic information beyond abundance alone and suggests that more diffuse immune infiltration may reflect more effective anti-tumor activity in ovarian cancer. Overall, FDA offers a flexible and interpretable framework for modeling spatial clustering across scales and identifying prognostic spatial features not captured by fixed-radius or distance analyses.

Dynamic models are central to the computational life sciences but typically contain unknown parameters that must be inferred from experimental data. High-throughput measurements have made this task increasingly challenging, yielding high-dimensional search spaces and non-convex objectives with many local optima. This makes the choice of optimisation method critical. However, existing empirical studies either consider only a limited number of benchmark problems or only a narrow spectrum of local, global and hybrid optimisation methods. Here, we present a comprehensive benchmark of a broad range of optimisation methods on a curated collection of parameter estimation problems, comprising 990 method-problem-pairs executed on two independent supercomputing infrastructures. Our evaluation quantifies success rates, solution quality and computational cost, revealing characteristic strengths and limitations of each approach. We find that optimisation methods separated into clear performance tiers. Building on these results, we implemented a new hybrid strategy that combines enhanced scatter search with the best-performing local solver, which showed robust performance and improved on the other scatter-search variants we tested. Our results provide practical guidance for selecting optimisation methods and thereby support more accurate and reliable model calibration.

Triple Negative Breast Cancer (TNBC) is the deadliest and most aggressive subtype of breast cancer, with poor prognosis and high rates of metastasis. Despite knowledge of metabolic rewiring in TNBC, the systems-level coordination of these adaptive pathways remains unmapped. This integrative systems-level analysis reveals key metabolic hub genes and identifies ATP1A2 as a significant prognostic marker. Analysis identified 764 differentially expressed genes, with 89 enriched biological processes predominantly involving metabolic pathways. Co-expression network analysis of 261 genes identified metabolic hub genes including LEP, ADIPOQ, and ATP1A2. To evaluate the prognostic framework, survival analysis of the top 10 hubs was performed on synthetic survival data, revealing ATP1A2 as a significant marker (p = 0.03) under Cox regression, with elevated expression associating with altered survival outcomes. By systematically mapping metabolic rewiring in TNBC, this work identifies ATP1A2 as an actionable therapeutic target and establishes a systems-level framework for rational drug discovery and patient stratification in this aggressive malignancy.

Cells adjust their internal circuits in response to changes in their environment. Hence, exposing cells to changing conditions provides a way to probe the intrinsic dynamics of cellular internal circuits. Metabolic networks are examples of such circuits since metabolic fluxes dynamically adjust when environmental conditions are transiently altered. Most existing theoretical frameworks focus on cellular metabolic steady states and do not consider the dynamics of changes in metabolic fluxes. In this work, we applied transfer function analysis from control theory to analyze the changes of NADH oxidative fluxes in the mitochondria and cytoplasm in mouse oocytes in response to dynamical perturbations of oxygen depletion and recovery. We observed an overshoot of NADH oxidative flux in the cytoplasm upon oxygen recovery which is absent in the mitochondrial NADH oxidative flux. Metabolic perturbation experiments and transfer function analysis indicate that this cytoplasmic NADH overshoot results from the coupling of the mitochondrial and cytoplasmic NADH cycles. The degree of overshoot is determined by competing timescales associated with the exchange rates of lactate and pyruvate with the media and their interconversion rates catalyzed by lactate dehydrogenase. Applying control theory to the data enables the inference of the exchange and conversion rates of pyruvate and lactate, allowing predictions of the contribution of lactate to mitochondrial respiration. Our work indicates that the oocytes maintain a homeostatic respiration rate across nutrient conditions by modulating the contribution of lactate to mitochondrial respiration.

Early-phase clinical trials of therapies for acute organ injury are typically small, uncontrolled, and must infer treatment activity using only tissue-damage biomarker changes over time. Interpretation is confounded by the temporal overlap of ongoing tissue damage and biomarker clearance. We address this with a Bayesian deconvolution framework that fits individual patient biomarker trajectories with an exponentially-modified Gaussian model. This model separates injury kinetics (peak release rate, injury duration, time to peak) from biomarker clearance, using a historic cohort as a null distribution. In a simulated phase 1 dose-escalation regenerative therapy trial, the framework reduces the minimum detectable treatment effect from 67.5% to 24.5% (a 2.76-fold improvement) and supports dose selection. Applied to published case-series data for another therapeutic, the framework recovers per-patient pharmacodynamic signatures consistent with pre-clinical mechanistic studies and supports smaller prospective clinical trial design. With indication-specific recalibration, the framework architecture conceptually transfers to other acute organ injuries where serum biomarkers reflect tissue damage. This offers a route to significantly reducing clinical trial sizes, supporting therapy dose-finding in non-controlled clinical trials, and identifying pharmacodynamic signatures.

Dynamic Flux Balance Analysis (DFBA) enables simulation of microbial culture dynamics under changing environmental conditions, but remains computationally expensive for tasks such as parameter calibration and fermentation optimization when applied using genome-scale metabolic models (GEMs). To address this challenge, we introduce Dynamic Flux Vector Balancing (DFVB), a reformulation of DFBA that solves an equivalent problem using a pre-computed, sparse basis of flux solutions that reduces the dimensionality of the internal optimization problem without information loss. Notably, DFVB provides a compact, interpretable representation of flux states that can readily identify dynamically inactive pathways and enable simulation-based automatic metabolic network reduction. We showed that DFVB produces the same culture dynamics as DFBA across multiple model scales and conditions, and identifies inactive reactions more accurately than Flux Variability Analysis (FVA) when compared to transcriptomic data profiles. Furthermore, computational performance analyses demonstrated that integrating DFVB with solver warm-start strategies and model reduction enhances computational efficiency relative to DFBA, yielding up to 3-fold reductions in simulation time for large-scale metabolic models. Finally, kinetic parameter estimation of culture dynamics with DFVB in two fermentation scenarios using a large-scale yeast GEM reached equal or higher prediction fidelity and narrower confidence intervals than DFBA, indicating improved parameter identifiability and robustness. Together, these results position DFVB as a scalable, robust, and biologically coherent framework for dynamic metabolic modeling, easing the integration of GEMs for culture dynamics simulation.

Type 2 diabetes is characterized by progressive aggregation of islet amyloid polypeptide (IAPP) within the islets of Langerhans, a process strongly implicated in beta-cell dysfunction and loss. Although oligomeric IAPP intermediates are widely considered the principal cytotoxic species, the relative contributions of the many biological and kinetic processes governing their formation, clearance, and conversion into fibrils remain poorly quantified. Here, a mathematical model of IAPP aggregation is developed that incorporates the physiology of beta-cell secretion and the microanatomy of the islet, including capillary-mediated clearance, enzymatic degradation, and the kinetics of oligomer and fibril formation within a well-mixed control volume. Building on the hypothesis that oligomers are the major cytotoxic species, the concept of accumulated cytotoxicity is introduced, defined as the time integral of the oligomer concentration, and a systematic sensitivity analysis of this quantity with respect to all model parameters is performed. The results reveal a striking hierarchy: only two parameters, the basal rate of IAPP monomer secretion and the rate constant for spontaneous oligomer dissociation, exert a first-order influence on long-term accumulated cytotoxicity, with dimensionless sensitivities approaching +1 and -1, respectively, while the effect of all other parameters remains subordinate and decays at long times. The model further shows that capillary clearance, owing to the physical exclusion of oligomers from fenestrated capillaries, selectively reduces fibril accumulation and amyloid deposition without affecting oligomer-mediated cytotoxicity, indicating that amyloid area fraction, the standard histological metric of disease severity, may not be a reliable surrogate for cytotoxic burden. The model predicts that approximately 48% of the islet area is replaced by amyloid after 30 years, broadly consistent with histological observations of advanced disease. These findings identify monomer secretion and oligomer dissociation as the most promising therapeutic targets to limit cytotoxic damage in type 2 diabetes and provide a quantitative framework for evaluating candidate intervention strategies.

Current approaches to cancer drug target discovery face two key limitations: poor translation of cell line-derived targets to patient tumors, and the lack of causal explanation of the regulatory mechanisms underlying target prioritization. Here we present BayesTx (Bayesian Therapeutics target discovery), a Bayesian network framework that integrates patient transcriptomics data with cell line data to identify causal therapeutic targets in cancer. BayesTx projects both data domains into a shared biological space of pathway and transcription factor activities, learns domain-specific causal graphs, and merges them through weighted edge aggregation with bootstrap consensus filtering. Do-simulation on the consensus network quantifies the causal effect of each transcription factor on cancer cell viability. Applied to breast cancer using TCGA-BRCA (The Cancer Genome Atlas breast cancer cohort) and DepMap (Cancer Dependency Map) datasets, the framework ranked 47 transcription factors by predicted causal impact, with gene-level targets further derived through regulon-based propagation. Top-ranked transcription factor (TF) targets were independently supported by survival analysis in external cohort data and pharmacogenomic drug response associations. Overall, BayesTx demonstrates that cross-domain Bayesian network modeling can bridge patient and cell line data to systematically identify causal therapeutic targets in cancer.

Cell states are increasingly conceptualized as attractors of high-dimensional dynamical systems, yet quantitative approaches for integrating phenotypic information into this framework remain limited. Here, we take an image-based approach that combines unsupervised machine learning (ML) with timelapse imaging to extract and characterize the temporal dynamics of morphological features. Using a cell line with endogenously tagged VE-cadherin, we acquired brightfield and fluorescence timelapse images of human induced pluripotent stem cell-derived endothelial cell (hiPSC-EC) monolayers, which adopt distinct phenotypes at two different magnitudes of shear stress in terms of their morphology, behavior, and VE-cadherin organization. To quantify these phenotypic cell states without segmentation, we trained a diffusion autoencoder to predict VE-cadherin signal from brightfield images. We identified interpretable ML-based features representing cell orientation, elongation, and local density. Treating these variables as dimensions of a morphological state space, we estimated a data-driven vector field and found that the two observed phenotypic cell states correspond to stable fixed points of the inferred dynamical system. Mapping measured cell migration coherence onto this space further distinguished the states. Imaging cells across intermediate shear stresses revealed a regime of bistability in which both states coexist, indicating that the shear-stress-dependent transition between endothelial cell states occurs as a bifurcation of the inferred dynamical system. Finally, we applied this method to study an N-terminal truncation of VE-cadherin, finding that mutated cells preserve alignment and coherent migration, but exhibit altered morphology and increased migration speed. This work demonstrates the applicability of a dynamical systems approach to quantitatively characterize morphological aspects of cell state from interpretable ML-based features.

Bulk RNA-seq remains the most accessible transcriptomic platform for tumor microenvironment (TME) characterization, yet existing computational approaches treat cell type abundance, pathway activity, and transcription factor (TF) activity as independent sources of information, missing the coordinated regulatory programs that define functional multicellular states. Here we introduce CellTFusion, a framework that integrates cell type deconvolution with transcriptional regulatory network analysis from bulk RNA-seq data to identify functional multicellular groups. CellTFusion produces a mixture representation of the TME in which each patient is described as a weighted combination of states, each capturing a distinct coordinated hallmark program. Applied to melanoma and bladder cancer cohorts, CellTFusion identified recurrent TME programs with opposing associations with immunotherapy responses that only emerged through joint multivariate modeling, and demonstrated superior cross-cohort transferability compared to established TME characterization tools.